# 相关库
from scipy import  stats
from statsmodels.graphics.tsaplots import plot_pacf as PACF   #偏自相关图
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.tsa.arima.model import ARIMA
import statsmodels.api as sm  # 统计相关的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import arch  # 条件异方差模型相关的库
import seaborn as sns     #seaborn画图
sns.set(color_codes=True) #seaborn设置背景

#①导入数据
data=pd.read_excel('中芯国际股价数据.xlsx')
data.set_index('date', inplace=True) #设定日期为索引
r2=np.log(data['close'])-np.log(data['close'].shift(1)) #计算对数收益率
r2=r2.dropna()  #去除包含NaN的行

#②使用偏自相关函数PACF对对数收益率序列定阶
#PACF图定阶
fig = PACF(r2, lags = 30)   
plt.show()

#信息准则AIC定阶
aicList = []
bicList = []
hqicList = []
for i in range(1,10):  #从1阶开始算
    od = (i,0,0)  # 这里使用了ARMA模型，order 代表了模型的(p,q)值，我们令q始终为0，就只考虑了AR情况。
    tempModel = ARIMA(r2,order=od).fit()
    aicList.append(tempModel.aic)
    bicList.append(tempModel.bic)
    hqicList.append(tempModel.hqic)
plt.figure(figsize=(15,6))
plt.plot(aicList,'r',label='aic value')
plt.plot(bicList,'b',label='bic value')
plt.plot(hqicList,'k',label='hqic value')
plt.legend(loc=0)

print(sm.tsa.arma_order_select_ic(r2,max_ar=6,max_ma=0,ic='aic')['aic_min_order'])
#返回（5，0），即AR模型定阶为5，通过PACF图也可以观察到第一个不在置信区间内的点横坐标为5

#③建立AR(5)模型
temp = np.array(r2) # 载入收益率序列
model =AutoReg(temp,lags=5)  
res = model.fit()  
out = 'AIC: {0:0.3f}, HQIC: {1:0.3f}, BIC: {2:0.3f}'
print(out.format(res.aic, res.hqic, res.bic))
print(res.summary())#输出模型结果

plt.rcParams['font.sans-serif'] = ['simhei'] #字体为黑体，正常显示中文字体
plt.rcParams['axes.unicode_minus'] = False #正常显示负号 

plt.figure(figsize=(10,4))
plt.plot(temp,'b',label='对数收益率') #核心画图代码 原始收益率序列
plt.plot(res.fittedvalues, 'r',label='AR model')#模型拟合结果
plt.legend()

#④模型检验：残差为白噪声则模型显著
#求残差序列
delta = res.fittedvalues  - temp[5:]  # 残差
plt.figure(figsize=(10,6))
plt.plot(delta,'r',label=' residual error')
plt.legend(loc=0)
plt.show()
acf,q,p = sm.tsa.acf(delta,nlags=12,qstat=True)  ## 计算自相关系数 及p-value
out = np.c_[range(1,13), acf[1:], q, p]
output=pd.DataFrame(out, columns=['lag', "AC", "Q", "P-value"])
output = output.set_index('lag')
print(output)
#p-value均大于0.05，在置信区间内，接受原假设：残差为白噪声

#⑤拟合优度
score = 1 - delta.var()/temp[5:].var()
print(score)
#很低，仅为0.01，效果不佳
